/*  CPLX.H - A complex structure definition and set of complex
    math operations for use with C.
    (C) 1998, Michael S Kelley   Last updated: 2-9-98

Data types:
    struct cplx - the complex data type, two members: re (real part),
                  and im (imaginary part)
Functions:
    struct cplx makecplx(double a, double b) - Returns a complex number a+bi

    double real(struct cplx z) - Returns the real part of z
    double imag(struct cplx z) - Returns the imaginary part of z

    struct cplx csum(struct cplx x, struct cplx y)  - Returns x+y
    struct cplx cdiff(struct cplx x, struct cplx y) - Returns x-y
    struct cplx cmult(struct cplx x, struct cplx y) - Returns x*y
    struct cplx cdiv(struct cplx x, struct cplx y)  - Returns x/y

    double modulus(struct cplx z) - Returns the modulus of z
    double phase(struct cplx z)   - Returns the angle of z (phi=arctan(imag/real))

    struct cplx cplxsqrt(struct cplx z)          - Returns the square root of z
    struct cplx conju(struct cplx z)              - Returns the conjuugate of z
    struct cplx cplxpow(struct cplx z, double k) - Returns z^k
    struct cplx cplxexp(struct cplx z)           - Returns e^z
*/

#ifndef __complex_mike
#define __complex_mike

#ifndef _MATH_H
#include <math.h>
#endif

struct cplx {
	double re;
	double im;
};

struct cplx makecplx(double a, double b);

double real(struct cplx z);
double imag(struct cplx z);

struct cplx csum(struct cplx x, struct cplx y);
struct cplx cdiff(struct cplx x, struct cplx y);
struct cplx cmult(struct cplx x, struct cplx y);
struct cplx cdiv(struct cplx x, struct cplx y);

double modulus(struct cplx z);
double phase(struct cplx z);

struct cplx cplxsqrt(struct cplx z);
struct cplx conju(struct cplx z);
struct cplx cplxpow(struct cplx z, double k);
struct cplx cplxexp(struct cplx z);

struct cplx makecplx(double a, double b)
{
	struct cplx result;
    result.re=a;
    result.im=b;
	return result;
}

double real(struct cplx z)
{
	return z.re;
}

double imag(struct cplx z)
{
	return z.im;
}

struct cplx csum(struct cplx x, struct cplx y)
{
	return makecplx(x.re+y.re, x.im+y.im);
}

struct cplx cdiff(struct cplx x, struct cplx y)
{
	return makecplx(x.re-y.re, x.im-y.im);
}

struct cplx cmult(struct cplx x, struct cplx y)
{
	return makecplx(x.re*y.re - x.im*y.im, x.re*y.im + x.im*y.re);
}

struct cplx cdiv(struct cplx x, struct cplx y)
{
	double d = y.re*y.re + y.im*y.im;
	return makecplx((x.re*y.re + x.im*y.im)/d, (-x.re*y.im + x.im*y.re)/d);
}

double modulus(struct cplx z)
{
	return sqrt(z.re*z.re + z.im*z.im);
}

double phase(struct cplx z)
{
	double d=atan(z.im/z.re);
	if (z.re < 0.0 && z.im <= 0.0) d=d+3.1415927;
	if (z.im > 0.0 && z.re <=  0.0) d=d+3.1415927; 
	return d;
	
}

struct cplx cplxsqrt(struct cplx z)
{
  double r, p;
  r = sqrt(modulus(z));
  p = phase(z)/2;

  return makecplx(r*cos(p), r*sin(p));
}

struct cplx conju(struct cplx z)
{
	return makecplx(z.re, -z.im);
}

struct cplx cplxpow(struct cplx z, double k)
{
  double r, p;
  r = pow(modulus(z), k);
  p = phase(z) * k;

  return makecplx(r*cos(p), r*sin(p));
}

struct cplx cplxexp(struct cplx z)
{
  return cmult(makecplx(exp(z.re), 0), makecplx(cos(z.im), sin(z.im)));
}  

#endif    /* __complex_mike  */
